/**
 *
 *                              Program/Module                                   
 *                                   from                                        
 *                 Selfverifying Solvers for Dense Systems of                    
 *                     Linear Equations Realized in C-XSC                        
 *                                                                               
 *                     Carlos Holbig and Walter Kraemer, 2003.                   
 *                                                                               
 *       Program developed in C-XSC by Bernardo Frederes Kramer Alcalde,         
 *           Paulo Sergio Morandi Junior and Carlos Amaral Holbig                
 *                                                                               
 * Purpose: Compute approximations to the solution x of Ax = b and an            
 *    approximations of the inverse of A. In both case the system (m x n) can be 
 *    square (m = n), over-determined (m > n) and under-determined (m < n).      
 */

#include "ManagerLss.h" 

/**
 * Construtor do objeto.
 * Parâmetros de entrada:
 *  rmatrix& A: matriz de coeficientes;
 *  rvector& b: matriz de termos independentes.
 * Parâmetros de saída:
 *  ivector& y: resultado intervalar retornado pelo solver.
 * 
 * Quando o solver for instanciado, as matrizes "user" e "sys" da classe são inicializadas
 * com as matrizes A, b e y passadas por parâmetro.
 * Além disso, o objeto parallelMachine contido na classe lss é instanciado.
 * 
 * Os membros userA, userB e userY são referências para os objetos de entrada do usuário. Os valores
 * desses objetos não são alterados no decorrer da execução do solver.
 * Os membros sysA, sysB e sysY são referências para as matrizes que o solver poderá alterar.
 * Qualquer processamento que é feito no solver, toma-se como base a referência dessas variáveis de sistema.
 */
ManagerLss::ManagerLss(rmatrix& A, rvector& b, ivector& y, int rank, int worldCommSize, int numThreads) : userA(A), sysA(A), userB(b), sysB(b), userY(y), sysY(y) {
	// Número de linhas e colunas do sistema
	this->nRows = Ub(this->userA,ROW) - Lb(this->userA,ROW) + 1;
	this->nCols = Ub(this->userA,COL) - Lb(this->userA,COL) + 1;
	
	this->setRank(rank);
	this->setWorldCommSize(worldCommSize);
	this->setNumThreads(numThreads);
}

/**
 * Destrutor do objeto.
 * 
 */
ManagerLss::~ManagerLss(){
}

/**
 * Retorna uma referência para a matriz de coeficientes do usuário.
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: rmatrix&
 * 
 */
inline rmatrix& ManagerLss::getA(void){
	return this->userA;
}

/**
 * Retorna uma referência para o vetor de termos independentes do usuário.
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: rvector&
 * 
 */
inline rvector& ManagerLss::getB(void){
	return this->userB;
}

/**
 * Retorna o vetor intervalar resultante do solver.
 * A referência retornada por esse método é a mesma retornada
 * pelo parâmetro de saída y do construtor do solver.
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: ivector&
 * 
 */
inline ivector& ManagerLss::getLssResult(void){
	return this->userY;
}

/**
 * Retorna o tipo de sistema que o usuário entrou.
 * Pode retornar os valores:
 * 	-1 = para sistemas sub-determinados (under-determined systems);
 * 	 0 = para sistemas quadrados (square systems);
 *	 1 = para sistemas sobre-determinados (over-determined systems).
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: int
 * 
 */
inline int ManagerLss::getSystemType(void){
	return this->systemType;
}

/**
 * Ajusta o tipo de sistema que o solver terá que tratar.
 * O sistema de entrada pode ser do tipo:
 * 	-1 = para sistemas sub-determinados (under-determined systems);
 * 	 0 = para sistemas quadrados (square systems);
 *	 1 = para sistemas sobre-determinados (over-determined systems).
 * 
 * Parâmetros de entrada:
 * 	int systemType, onde systemType = [-1|0|1]
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: void
 * 
 */
void ManagerLss::setSystemType(int systemType){
	if(systemType < -1 || systemType > 1){
		this->systemType = 0;
	} else {
		this->systemType = systemType;
	}
}

/**
 * Dispara a execução do solver.
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: int
 * 
 */
int ManagerLss::startLss(void){
	int errcode = 0;

	if(nRows != Ub(this->userB) - Lb(this->userB) + 1)
		return(3);  // b : wrong dimension

	if(nCols != Ub(this->userY) - Lb(this->userY) + 1)
		return(4);  // y : wrong dimension

	if(nRows == nCols){
		this->setSystemType(0);
	} else {
		if(nRows > nCols){
			this->setSystemType(1);
			this->over_lss(); // over-determined system
		} else {
			this->setSystemType(-1);
			this->under_lss(); // under-determined system
		}
	}
	
	errcode = this->square_lss(); // square system

	return errcode;
}

/**
 * Executa o solver para o sistema linear quadrado.
 * O resultado contido em sysY é uma aproximação da solução do sistema Ax = b 
 *
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	nenhum.
 * Valor de retorno: int, onde pode assumir os valores:
 * 	0: sysY é a aproximação da solução;
 * 	1: nenhuma aproximação obtida (bad condition);
 * 	2: matrix A é singular (aproximação não obtida).
 */
int ManagerLss::square_lss(void){
	// Código de retorno do método
	int errcode = 0;
	// Dimensão do sistema
	int dim = Ub(this->sysA,ROW) - Lb(this->sysA,ROW) + 1;
	
	// Treat trivial case separately
	if(dim == 1) {
		idotprecision Accu;

		if(this->sysA[Lb(this->sysA, ROW)][Lb(this->sysA, COL)] == 0.0){
			errcode = 2;
		} else {
			this->sysY[Lb(this->sysY)] = interval( this->sysB[Lb(this->sysB)] ) / this->sysA[Lb(this->sysA,ROW)][Lb(this->sysA,COL)];
			errcode = 0;
		}
	
	// when dim != 1
	} else {
		rmatrix R1(Lb(this->sysA,ROW), Ub(this->sysA,ROW), Lb(this->sysA,ROW), Ub(this->sysA,ROW));
		rmatrix R2(Lb(this->sysA,ROW), Ub(this->sysA,ROW), Lb(this->sysA,ROW), Ub(this->sysA,ROW));
		
		// Chama o método para o cálculo em precisão simples
		errcode = this->compute_single_r(R1, R2);

		// Se o resultado não foi obtido, chama o método para o cálculo em dupla precisão
		if (errcode == -1) {
			errcode = this->compute_double_r(R1, R2);
		}
	}
	
	// Trata os casos em que o sistema não é quadrado
	switch(this->getSystemType()){
		case -1:// Under-determined system
			this->userY = this->sysY(this->nRows+1, this->nRows + this->nCols);
			break;
		case 1: // Over-determined system
		 	this->userY = this->sysY(1, this->nCols);
		 	break;
	}
	
	return errcode;
} // end SQUARE_LSS


/**
 * Linear system: over-determined case
 * The result Y is an enclosure of the solution x of AH*A*x = AH*b,
 * i.e. x is least squares solution of Ax = b.
 *                                                 | A  -I | |x|   |b|
 * From Ax = b we generate the (n+m)x(n+m)-system  |       | | | = | |
 *                                                 | 0  AH | |y|   |0|
 * 
 * Here, AH is the hermitian of A (transpose in the real case)
 * 
 */
void ManagerLss::over_lss(void) {
	int i = 0, j = 0;

	int dim = this->nRows + this->nCols;

	rmatrix BIG_A(1,dim,1,dim);
	rvector BIG_b(1,dim);
	ivector BIG_y(1,dim);

	BIG_A(1, this->nRows, 1, this->nCols) = this->userA;

	// BIG_A( 1,nRows, nCols+1,nCols+nRows ) = -Id(nRows);
	for(i = 1; i <= this->nRows; i++){
		for(j = this->nCols + 1; j <= this->nCols + this->nRows; j++){
			(j == i+2)? BIG_A[i][j]=-1 : BIG_A[i][j] = 0;
		}
	}

	BIG_A(this->nRows+1, this->nRows + this->nCols, 1, this->nCols) = 0.0;
	BIG_A(this->nRows+1, this->nRows + this->nCols, this->nCols + 1, this->nCols + this->nRows) = transp(this->userA);
	BIG_b(1, this->nRows) = this->userB;
	BIG_b(this->nRows+1, this->nRows + this->nCols) = 0.0;
	
	this->sysA = BIG_A;
	this->sysB = BIG_b;
	this->sysY = BIG_y;

} // end OVER_LSS

/*
 * Linear system: under-determined case
 * The result Y is an enclosure of Y = AH*x with A*AH*x = b,
 * i.e. y is solution of Ay = b with minimal Euklidian norm.
 *
 *                                                 |AH  -I | |x|   |0|
 * From Ax = b we generate the (n+m)x(n+m)-system  |       | | | = | |
 *                                                 | 0   A | |y|   |b|
 * 
 * Here, AH is the hermitian of A (transpose in the real case)
 * 
 */
void ManagerLss::under_lss(void) {
	int i = 0, j = 0;

	int dim = this->nRows + this->nCols;

	rmatrix BIG_A(1,dim,1,dim);
	rvector BIG_b(1,dim);
	ivector BIG_y(1,dim);

	BIG_A(1, this->nCols, 1, this->nRows) = transp(this->userA);

	//BIG_A( 1,n, m+1,m+n ) = -Id(m);
	for(i = 1; i <= this->nCols; i++){
		for(j = this->nRows + 1; j <= this->nCols + this->nRows; j++){
			(j == i+2)? BIG_A[i][j]=-1 : BIG_A[i][j] = 0;
		}
	}

	BIG_A(this->nCols + 1, this->nCols + this->nRows, 1, this->nRows) = 0.0;
	BIG_A(this->nCols + 1, this->nCols + this->nRows, nRows + 1, this->nRows + this->nCols) = this->userA;
	BIG_b(1, this->nCols) = 0.0;
	BIG_b(this->nCols + 1, this->nCols + this->nRows) = this->userB;

	this->sysA = BIG_A;
	this->sysB = BIG_b;
	this->sysY = BIG_y;
} // end UNDER_LSS

/**
 * Calcula uma aproximação em precisão simples da inversa de A.
 * Utiliza como base os sistemas referenciados por sysA e sysB.
 * 
 * Parâmetros de entrada:
 * 	nenhum.
 * Parâmetros de saída:
 * 	rmatrix& R1
 * 	imatrix& C
 * Valor de retorno: int, onde pode assumir os valores:
 * 	0: sysY é a aproximação da solução;
 * 	1: nenhuma aproximação obtida (bad condition);
 * 	2: matrix A é singular (aproximação não obtida).
 * 
 */
int ManagerLss::compute_single_r(rmatrix& R1, rmatrix& R2){
	rvector x0(Lb(this->sysA,1),Ub(this->sysA,1)), 
			x1(Lb(this->sysA,1),Ub(this->sysA,1)),
			y0(Lb(this->sysA,1),Ub(this->sysA,1)),
			x1temp(Lb(this->sysA,1),Ub(this->sysA,1));

	ivector Y1(Lb(this->sysA,1),Ub(this->sysA,1)),
			YA(Lb(this->sysA,1),Ub(this->sysA,1)),
			Z(Lb(this->sysA,1),Ub(this->sysA,1)),
			Y1temp(Lb(this->sysA,1),Ub(this->sysA,1));
	
	dotprecision accu;
	idotprecision iaccu;
	
	const real sqrt_01 = 0.31622777; // raiz quadrada de 0,1, but why square root of 0,1?????
	const real delta = 1e-15;
	
	int errcode = -1;
	int err = 0;
	
	cout << "Iniciando inversao de A" << endl;
	// Compute an approximate inverse R of sysA using the Gauss-Jordan Algorithm
	gaussJordan(this->sysA, R1);

	// floating point defect iteration: result is x1
	if (err == 0) {
		// Compute an approximation x' = R * sysB of x
		x1 = R1 * this->sysB;

		real bound = 100.0 * sqrt_01;
		real p = 0.0;
		int k = 0;

		// improve x' by an iterated defect correction
		do { // iterate x = x + R*(b-Ax)
			k = k + 1;

			x0 = x1; // x1 contém o resultado parcial do sistema

			// x1 := #*(b - A*x0)
			// Determina os valores do erro com relação à b (b - A * x), pois A*x = b
			for (int i = Lb(x0); i <= Ub(x0); i++) {
				// Percorre o vetor sysB do sistema jogando cada elemento para um acumulador dotprecision
				accu = this->sysB[i];
				// Determina o erro do resultado, com relação ao vetor de termos independentes sysB
				// ou seja, multiplica-se A e x (onde x é a aproximação do resultado) e subtrai-se o resultado da multiplicação
				// do respectivo valor de b. O vetor x é uma aproximação, logo o resultado de A*x também é aproximado.
				// O erro dessa aproximação é determinado subtraindo-se do valor correto (contido em b) o valor aproximado.
				accumulate(accu, -this->sysA[Row(i)], x0); // accu -= this->sysA[Row(i)] * x0;
				
				// Ao final do laço, x1 contém os valores de erro de cada posição do vetor
				x1[i] = rnd(accu);
			}
			
			// x1 contém os valores do erro absoluto com relação à b
			// x1 := #*(x0 + R1*x1)
			for (int i = Lb(x1); i <= Ub(x1); i++) {
				// Para cada valor do resultado aproximado contido em x0, joga-o para o acumulador dotprecision
				accu = x0[i];			
				// Realiza o produto escalar ótimo de cada linha da inversa aproximada R1 de sysA
				// pelo vetor de erros x1, determinado no laço anterior
				// O produto escalar do vetor que contém o erro absoluto com a linha i da matriz inversa
				// resulta no escalar que deve ser somado ou subtraído do resultado aproximado.
				accumulate(accu, R1[Row(i)], x1); // accu += linha de R1 * valores de erro
				
				// Ao final do laço, x1temp contém os valores do resultado aproximado somados do produto escalar ótimo
				// do valor do erro pela inversa aproximada.
				// Ou seja, x1temp vai ter resultado aproximado + linha da inversa * erro do resultado aproximado com relação à b
				x1temp[i] = rnd(accu);
			}
			
			// Joga para x1 os valores do vetor x1temp determinado no último laço
			x1 = x1temp;
			
			// Determina o maior erro relativo de x1 com relação à x0 e retorna para p
			p = rel(x1, x0);
			
			if(k > 5) {
				bound = bound * sqrt_01;
			}
		
		// Iterar enquanto o maior erro relativo (p) for menor que 31.622777 (bound) ou
		// k for menor ou igual à 5. E enquanto o maior erro relativo for maior que 0.000000000000001 (delta).	
		} while ((p < bound || k <= 5) && p >= delta);
		
		// x0 é alterado
		guess_zeroes(x1, x0);
		
		// Compute enclosure y0 + y1 of the residuum b-sysA*x1 of the aproximation x1
		// and initialize Y1:=Z:= R1*(b-sysA*x1), C:= I-R1*sysA
		
		// y0 := #*(b-sysA*x1) and Y1 := ##(b-sysA*x1-y0)
		for (int i=Lb(x1); i <= Ub(x1); i++) {
			accu = this->sysB[i];
			accumulate(accu, -this->sysA[Row(i)], x1);
			rnd(accu, y0[i]);
			
			accu = accu - y0[i];
			rnd(accu, Y1[i]);
		}
			
		// Y1 := ##( R1*y0 + R1*Y1 );
		for (int i=Lb(R1,1); i <= Ub(R1,1); i++) {
			iaccu = 0.0;
			accumulate(iaccu, R1[Row(i)], y0);
			accumulate(iaccu, R1[Row(i)], Y1);
			rnd(iaccu, Y1temp[i]);
		}
			
		Y1 = Y1temp;
		Z = Y1;
		
		if (Z == null(Z)) {
			this->sysY = x1; // exact solution! (however, not necessarily unique!)
			errcode = 0;
		} else {
			imatrix C(Lb(this->sysA,ROW), Ub(this->sysA,ROW), Lb(this->sysA,ROW), Ub(this->sysA,ROW));
			
			cout << "Calculando IminusAB paralelo" << endl;
			IminusAB(R1, this->sysA, C, R2, this->getRank(), this->getNumThreads(), this->getWorldCommSize());
			
			// interval iteration until inclusion is obtained (or max. iteration count)
			k = 0;
			real eps1 = 0.1;
			bool ready = false;
			
			do {
				if (k >= 5){
					eps1 = 5 * eps1;
				}
				
				k = k + 1;
				YA = Blow(Y1, eps1);
				Y1 = Z + C * YA;
				ready = in(Y1, YA); // Verifica se Y1 está em YA
			} while ( !ready && k < 10 && !too_bad(Y1) );

			// output of the result
			if (ready) {
				this->sysY = x1 + Y1;
				errcode = 0;
			}
		}
	} // end USE_SINGLE_R
	
	return errcode;
} // end lss::compute_single_r

/**
 * Calcula uma aproximação em precisão dupla da inversa de A.
 * Utiliza como base os sistemas referenciados por sysA e sysB.
 * 
 * Parâmetros de entrada:
 * 	R1: inversa de A.
 * Parâmetros de saída:
 *  nenhum.
 * Valor de retorno: int, onde pode assumir os valores:
 * 	0: sysY é a aproximação da solução;
 * 	1: nenhuma aproximação obtida (bad condition);
 * 	2: matrix A é singular (aproximação não obtida).
 * 
 */
int ManagerLss::compute_double_r(rmatrix& R1, rmatrix& R2){
	rmatrix D(Lb(this->sysA,1),Ub(this->sysA,1),Lb(this->sysA,2),Ub(this->sysA,2)), 
			R2temp(Lb(this->sysA,1),Ub(this->sysA,1),Lb(this->sysA,2),Ub(this->sysA,2)),
			R2aux(Lb(this->sysA,1),Ub(this->sysA,1),Lb(this->sysA,2),Ub(this->sysA,2));

	rvector x0(Lb(this->sysA,1),Ub(this->sysA,1)), 
			x1(Lb(this->sysA,1),Ub(this->sysA,1)),
			y0(Lb(this->sysA,1),Ub(this->sysA,1)),
			x1temp(Lb(this->sysA,1),Ub(this->sysA,1));

	ivector Y1(Lb(this->sysA,1),Ub(this->sysA,1)),
			YA(Lb(this->sysA,1),Ub(this->sysA,1)),
			Z(Lb(this->sysA,1),Ub(this->sysA,1)),
			Y1temp(Lb(this->sysA,1),Ub(this->sysA,1));
	
	const real sqrt_01 = 0.31622777;
	const real delta = 1e-15;

	int k = 0, err = 0, errcode = -1;

	dotprecision accu;
	idotprecision iaccu;

	// R1 contém a inversa de sysA obtida em compute_single_r
	// R2 conterá o resultado da multiplicação de sysA com sua inversa
	// O resultado da multiplicação de uma matriz por sua inversa resulta na matriz identidade:
	
	cout << "Iniciando inversao de R2" << endl;
	// Calcula a inversa da matriz identidade obtida
	gaussJordan(R2, R2);
	
	if (err == 0) { // floating point defect iteration: result is x1+x0
		cout << "Calculando DR2temp paralelo" << endl;
		calculateDR2temp(R2, R1, R2temp, D, this->getRank(), this->getNumThreads(), this->getWorldCommSize());

		R2 = R2temp;
		R1 = D;
		
		real bound = 100.0 * sqrt_01;
		real p = 0.0;

		// x1 := #*( R1*b + R2*b ) e x0 = #*( R1*b + R2*b - x1):
		for (int i=Lb(R1,1);i<=Ub(R1,1);i++) {
			accu = 0.0;
			accumulate(accu,R1[Row(i)],this->sysB);
			accumulate(accu,R2[Row(i)],this->sysB);
			rnd(accu, x1[i]);
			
			accu = accu - x1[i];
			rnd(accu, x0[i]);
		}
		
		k = 0;
		
		do { // iteration x = x + (R1+R2)*(b-Ax), x = x1 + x0)
			k = k+1;
		    
			// y0 = #*(b - A*x1 - A*x0)
			for (int i=Lb(this->sysA,1);i<=Ub(this->sysA,1);i++) {
				accu = this->sysB[i];
				accumulate(accu,-this->sysA[Row(i)],x1);
				accumulate(accu,-this->sysA[Row(i)],x0);
				y0[i] = rnd(accu);
			}
			
			// y0 := #*(x0 + R1*y0 + R2*y0);
			for (int i=Lb(R1,1);i<=Ub(R1,1);i++) {
				accu = x0[i];
				accumulate(accu,R1[Row(i)],y0);
				accumulate(accu,R2[Row(i)],y0);
				x1temp[i] = rnd(accu);
			}

			y0 = x1temp;
		    
			p = rel(x1+y0,x1+x0);
			y0 = x1 + y0;
		    
			// x0 := #*(x1 + x0 - y0)
			for (int i=Lb(x1);i<=Lb(x1);i++) {
				accu = x1[i] + x0[i] - y0[i];
				x0[i] = rnd(accu);
			}
			
			x1 = y0;
		    
			if (k>5){
				bound = bound * sqrt_01;
			}
		} while ( (p<bound || k<=5) && p>=delta );
		
		// compute enclosure y0+Y1 of the residuum b-A*x1 of the approximation x1 and
		// initialize Y1:= Z:= (R1+R2)*(b-A*x1), C:= I-(R1+R2)*A

		// y0 := #*(b-A*x1) e Y1 := ##(b-A*x1-y0)
		for (int i=Lb(x1);i<=Ub(x1);i++) {
			accu = this->sysB[i];
			accumulate(accu,-this->sysA[Row(i)],x1);
			rnd(accu, y0[i]);
			
			accu = accu - y0[i];
			rnd(accu, Y1[i]);
		}

		// Y1 := ##(R1*y0 + R2*y0 + R1*Y1 + R2*Y1 );
		for (int i=Lb(R1,1);i<=Ub(R1,1);i++) {
			accu = 0.0;
			accumulate(accu,R1[Row(i)],y0);
			accumulate(accu,R2[Row(i)],y0);
			iaccu = accu;
			accumulate(iaccu,R1[Row(i)],Y1);
			accumulate(iaccu,R2[Row(i)],Y1);
			rnd(iaccu,Y1temp[i]);
		}
			
		Y1 = Y1temp;
		Z = Y1;
		
		if (Z == null(Z)) {
			this->sysY = x1; // exact solution! (however, not necessarily unique!)
			errcode = 0;
		} else {
			imatrix C(Lb(this->sysA,ROW), Ub(this->sysA,ROW), Lb(this->sysA,ROW), Ub(this->sysA,ROW));

			cout << "Calculando IminusABminusCB paralelo" << endl;
			IminusABminusCB(R1, R2, this->sysA, C, this->getRank(), this->getNumThreads(), this->getWorldCommSize());
			
			// interval iteration until inclusion is obtained (or max. iteration count)
			k  = 0;
			real eps1 = 0.1;
			bool ready = false;

			do {
				if (k >= 5){
					eps1 = 5*eps1;
				}
				k = k + 1;
				YA = Blow(Y1, eps1);
				Y1 = Z + C * YA;
				ready = in(Y1,YA);
			} while ( !ready && k < 10 && !too_bad(Y1) );

			// output of the result
			if (ready) {
				this->sysY = x1 + Y1;
				errcode = 0;
			} else {
				errcode = 1;
			}
		}
	} else {
		errcode = 2;
	}
		
	return errcode;
} // end lss::compute_double_r

int ManagerLss::getRank(void){
	return this->rank;
}

void ManagerLss::setRank(int rank){
	if(rank > -1){
		this->rank = rank;
	} else {
		this->rank = 0;
	}
	
	return;
}

int ManagerLss::getWorldCommSize(void){
	return this->worldCommSize;
}

void ManagerLss::setWorldCommSize(int worldCommSize){
	if(worldCommSize > 1){
		this->worldCommSize = worldCommSize;
	} else {
		this->worldCommSize = 1;
	}
	
	return;
}

int ManagerLss::getNumThreads(void){
	return this->numThreads;
}

void ManagerLss::setNumThreads(int numThreads){
	if(numThreads > 1){
		this->numThreads = numThreads;
	} else {
		this->numThreads = 1;
	}
	
	return;
}
